Spatial Analysis and Statistical Modeling with R and spmodel - 2
2024 Society for Freshwater Science Conference
Ryan Hill
EPA (USA)
Michael Dumelle
EPA (USA)
2024-06-02
GIS in R
Maintaining all analyses within a single software (R) can greatly simplify your research workflow. In this section, we’ll cover the basics of doing GIS in R.
Goals and Motivation
Understand the main features and types of vector data.
Generate point data from a set of latitudes and longitudes, such as from fields sites.
Read, write, query, and manipulate vector data using the sf package.
r <-rast(ncol=10, nrow =10)r[] <-runif(n=ncell(r))r
class : SpatRaster
dimensions : 10, 10, 1 (nrow, ncol, nlyr)
resolution : 36, 18 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
source(s) : memory
name : lyr.1
min value : 0.01459834
max value : 0.99461202
Raster data
plot(r)
Basic raster in R.
Raster data
We can access data from specific locations within a raster
# Access data from the ith location in a rasterr[12]
lyr.1
1 0.4480738
# Access data based on row and columnr[2, 2]
lyr.1
1 0.4480738
Raster data
Rasters can be stacked which can make them very efficient to work with
# Create 2 new rasters based on raster rr2 <- r *50r3 <-sqrt(r *5)# Stack rasters and rename to be uniques <-c(r, r2, r3)names(s) <-c('r1', 'r2', 'r3')
Raster data
plot(s)
Working with real data
There are several packages for accessing geospatial data from the web
We will use the FedData package, but numerous other packages exist to access data within and without the U.S.
One useful example is the elevatr package for accessing elevation data around the world
Working with real data
We will walk through an example of extracting information from a raster using a polygon layer. To do this we will:
Select just Corvallis among oregon_cities
Add a 10,000m buffer
Download National Elevation Data
Transform the projection system of the elevation data to match Corvallis
Calculate the average elevation within 10km of Corvallis
Working with real data
Buffer Corvallis
library(FedData)# Select just Corvallis and calculate a 10,000-m buffercorvallis <- oregon_cities %>%filter(cities =='Corvallis') %>%st_buffer(10000)
Working with real data
Download NED based on Corvallis buffer
# Download national elevation data (ned)ned <- FedData::get_ned(template = corvallis,label ="corvallis")
Transform the CRS
ned <- terra::project(ned, 'epsg:5070',method ='bilinear')
Working with real data
Mean elevation within polygon
# zonal function in terra to calculate zonal statisticsterra::zonal(ned, # Need to convert corvallis `sf` object to terra vector terra::vect(corvallis), # Metric to be calculated mean, na.rm = T)
Layer_1
1 119.6678
Your Turn
Read in U.S. cities with data('us.cities') from the maps library
Select the city of your choice and buffer it by 10Km. (We suggest converting to an equal area projection first)
Read in National Elevation Data for your city with the FedData package
Transform CRS of elevation data to match city
Calculate the mean elevation within 10km of your city
Characterizing watersheds is fundamental to much of our work in freshwater science
Although it is more than we can cover in today’s workshop, we want you to be aware that there are several options for delineating watersheds in R
We’ll provide two examples of how to delineate watersheds within the conterminous U.S. using two online services
USGS StreamStats
The USGS’s StreamStats is an online service and map interface that allows users to navigate to a desired location and delineate a watershed boundary with the click of a mouse:
nhdplusTools is an R package that can access the Network Linked Data Index (NLDI) service, which provides navigation and extraction of NHDPlus data: https://doi-usgs.github.io/nhdplusTools/
nhdplusTools includes network navigation options as well as watershed delineation
The delineation method differs from StreamStats (based on National Hydrography IDs)
nhdplusTools
library(nhdplusTools)# Simple feature option to generate point without any other attributescal_pt <-st_sfc(st_point(c(longitude, latitude)), crs =4269)# Identify the network location (NHDPlus common ID or COMID)start_comid <- nhdplusTools::discover_nhdplus_id(cal_pt)# Combine info into list (required by NLDI basin function)ws_source <-list(featureSource ="comid", featureID = start_comid)cal_ws2 <- nhdplusTools::get_nldi_basin(nldi_feature = ws_source)